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Abstract Numerical computation of the time evolution of the mass transfer rate in a close binary can be and, in particular, 
has been a computational challenge. Using a simple physical model to calculate the mass transfer rate, we show that for a 
simple explicit iteration scheme the mass transfer rate is numerically unstable unless the time steps are sufficiently small. In 
general, more sophisticated explicit algorithms do not provide any significant improvement since this instability is a direct 
result of time discretization. For a typical binary evolution, computation of the mass transfer rate as a smooth function of time 
limits the maximum tolerable time step and thereby sets the minimum total computational effort required for an evolutionary 
computation. By methods of "Controlling Chaos" it can be shown that a specific implicit iteration scheme, based on Newton's 
method, is the most promising solution for the problem. 
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1. Introduction 

To compute the long-term evolution of a semi-detached binary, 
the numerically determined mass transfer rate, which is a func- 
tion of current stellar and binary parameters, must be a suffi- 
ciently smooth function of time. The simplest approach to ob- 
tain the mass M 2 of the mass losing secondary star is an explicit 
forward integration in time: 

M 2 (f n+ i) = M 2 (t n ) + M 2 (t n ) At. (1) 

Here, Af = t n+ \ - t n denotes the length, and f n and t n+ i the 
start and end points of the nth time step. Every change of M 2 
within one time step causes changes of other stellar and binary 
parameters which affect M 2 , and an additional change of M 2 
would be necessary to obtain a consistent result for the current 
time step. This feedback between M 2 and other stellar and bi- 
nary parameters, especially the radius R 2 and the critical Roche 
radius Rr 2 of the donor star, can destabilize the numerically 
computed mass transfer rate: M 2 can become non-continuous 
between successive time steps and can even show strong fluc- 
tuations around its secular mean value. These numerical effects 
have been known for quite a long time (numerical experience 
by the authors; U. Kolb, K. Schenker, priv. comm.). Not sur- 
prisingly, only very few evolutionary calculations that show 
these instabili ties of the m ass transfe r rate have be en published, 
e.g..lKolb & RitteiUl99(i Fig. 3-5>.ISarnal dl992i Fig. 1 & 9), 
D'AntonaUl994lFig . 2). iKolb et all d2000l Fig. 3 & 4), and 
Schenker et all (120021 Fig. 5). 
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Attempts to use a different explicit integration scheme that 
takes into account not only M 2 (t n ) but also M 2 (t n -j) for certain 
values of j (e.g., for a particular average over the last m time 
steps) did not show major improvements. The reason for this 
behaviour of the mass transfer rate has remained unknown. 

The main purpose of this paper is to show analytically why 
these numerical instabilities exist, what they are, which meth- 
ods are suitable to suppress them, and which ones are not. Even 
in the case of the proposed implicit algorithm, calculation of 
the mass ransfer rate still limits the maximum tolerable time 
step in a numerical computation and thus sets the minimum to- 
tal computational effort required for carrying out such a binary 
evolution. To illustrate this point: by imposing a fixed mass loss 
rate on a single low-mass main sequence star, up to 10 % of the 
total mass of the star can be removed per time step (eventually 
after 1 or 2 initial time steps with a lower rate), and the stel- 
lar model still converges. But if the mass loss rate is coupled 
to the binary parameters, even in the case of our proposed im- 
plicit iteration scheme, typically no more than a few 10~ 3 of 
the stellar mass can be removed per time step without losing 
convergence 1 . When using the explicit algorithm Q, the cor- 
responding value is much lower: often only about 10~ 5 - 10 
or even less of the total mass can be removed per time step if 
fluctuations of the mass transfer rate by up to several orders of 
magnitude are to be avoided. 



1 This upper limit of about 10~ 3 varies and depends on the starting 
value for the iteration, on the "smoothness" of the stellar input physics, 
and on the "smoothness" of the stellar structure. In principle, this is not 
a limitation of the implicit mass transfer algorithm. 
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This paper is organized as follows: First, we discuss in 
Sect. |2]the necessary input physics before we show in Sect. [5] 
that the mass transfer rate in our model is physically stable. 
Subsequently, in Sect. 14.11 we prove mathematically that for 
a simple explicit algorithm the time-discretized mass transfer 
rate becomes unstable if At is greater than a critical value. 
We also give an estimate of how many time steps are re- 
quired to calculate the binary evolution. Next, we show in 
Sect. l4.2l that for increasi ng At the mass tran s fer rat e undergoes 
a Feigenbaum scenario (iFei ^en b aumlll97 8i Il98(l for a more 
accessible review see, e.g.. lThompson & Stewarl 19861) and fi- 
nally becomes chaotic. In Sect. l4.3l we discuss how the onset of 
chaos can be suppressed in terms of "Controlling Chaos" and 
we present an implicit iteration scheme for the integration of 
the mass transfer rate. Finally, in Sect. [5] we discuss a number 
of points which have to be considered for a practical implemen- 
tation. 



2. Input physics 

In this paper we u se the same nomenclature as in 
Buning & Ritterl J2004 . hereafter called Paper I. For the mass 
transfer rate we use 



(AR\ 

M 2 = -Mo expl — I, 



where 



AR '.— R2 — Rr,2 



(2) 



(3) 



denotes the difference between the radius Ri and the Roche 
radius R^2 of the donor, Hp the photospheric pressure scale 
height, and Mq > a weakl y varying fu nction of several system 
parameters (for details, see lRitterll9 88). 

Instead of (0 other relations betwee n M? and AR hav e also 
been used in the literature. For example iTout et alJdl997h have 
adopted a power-law dependence M2 x (AR/R2) i . However, 
a mass transfer prescription other than results mainly in a 
different characteristic scale length 



H :-- 



/ dln(-M 2 ) 
dAR 



(4) 



at the secular mass transfer rate M2 which is given by 



3. The time-continuous system 

The mass transfer rate M2 is uniquely invertibly coupled to AR 
by @. For simplicity, we will consider the time evolution of 
AR instead of M2 which is given by Eq. [36] of Paper I: 



±M = (ft 



M 2 (AR) R 2 



(6) 



Since for our stability considerations we neither take into ac- 
count changes in the system parameters (£ s , £r, tV) nor irradi- 
ation of the donor as in Paper I, our model is completely de- 
scribed by this 1 -dim. autonomous differential equation which 
simplifies the linear stability analysis significantly. 

At the stationary value AR which is the only fixed point 
(FP) of (jS) and which is equivalent to the stationary mass trans- 
fer rate, using (0 we get 



— dF(AR) 

DF(AR) := — - 

V ; dAR 



(f.-fn) 



R 2 M 2 

Th~M~2 



(7) 



(cf. Eq. [44] of PaperJ). Since M 2 < and ^ s > £ R , DF is 
negative, not only for AR, but even for all AR E R. This means 
that the stationary mass transfer rate, i.e., AR is stable and that 
all solutions AR(t) of converge 2 to AR. The convergence 
occurs on a time scale of 



Hp 
R 2 



AR{t) = AR ± [AR(0) - AR] e 



DF t 



(8) 



(9) 



in the lin earized system. This h as also been discussed in more 
detail bv lD'Antona et alJ dl989h . 



4. The time-discretized system 

4.1. The explicit algorithm 

The simplest way to obtain the donor mass numerically as a 
function of time is an explicit forward integration of M2 as 
given by Q. In this time-discretized system the evolution of 
AR is given by an iteration equation which follows directly 
from 



A/? 



n+l 



AR n 



R2 



ft) ^-M 2 (AR n ) 



At 



(10) 



Mt_ 

H 2 



(5) 



Here, tm is the secular time scale of the mass loss, t' a the driv- 
ing time scale including thermal relaxation of the donor, £ s and 
£r are respectively the adiabatic mass radius exponent of the 
donor and the mass radius exponent of the Roche radius, where 
£ s > £r is required for physical stability of mass transfer (for 
details, see Paper I). For the mass transfer prescription ©, we 
have H = Hp. 



Here, AR n denotes the value of AR at t — t n which is in practice 
a (numerical) approximation for the "real" value of AR in the 
time-continuous system. 

Obviously, (|6j and dl0> . i.e., the time-continuous and the 
time-discretized systems have the same FP since F(AR) — if 
and only if Q>(AR) = AR. But, since <&(AR) is a map while 



2 This can be proven by using the fact that F is continuous and has 
only one FP. Thus, DF(AR)_<_0 for all AR e R together with F(AR) = 
implies F % <=> AR S AR so that from J6j the convergence of all 
solutions to AR can be concluded. 
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F(AR) is a vector field 3 , the condition for stability is dif- 
ferent: according to the Hartm ann-Grobmann theorem (e.g., 
iGuckenheimer & Holmesiri983l) . a FP x of a map O is stable 
if the absolute values of all eigenvalues of the Jacobi matrix 
DO of O at x are less than unity, and x is unstable if at least 
one eigenvalue has an absolute value greater than unity. 

In the case of <f>(AR) from JlOi this means: AR is stable if 



|D(D| := 



dO 




dAR 





R 2 M 2 
1 + U s - £r) — — -At 
h h Hp M 2 



< 1, 



(11) 



and unstable if the absolute value is greater than unity. As can 
be easily shown by using Q, Eq. il It is equivalent to 



A? < 2— t' 

R 2 d 



At n 



(12) 



This means that in the neighbourhood of AR all solutions of the 
explicit iteration scheme converge to the FP if Af is less than 
the critical time step length Af max . But for Af > Af max the FP is 
unstable and the solutions of the time-discretized system i ll Oi l 
diverge although the solutions of the time-continuous system 
(j6jl converge to the FP. 

If the initial value of the iteration is ARq, then in the lin- 
earized system the orbit of A/?n, i.e., {A/?o, AR\, AR 2 , . . .} is 
given by 

AR n+ i = <D(AR n ) ~AR + [D^(AR)f +l (AR -AR\, (13) 

which can be shown by induction over n. 

The following cases are possible (without proof): 

1. Af < |f t^: the orbit converges directly to the FP; for Af — > 
the continuous limit is reached. 

2. At - ^t' a : the orbit converges after one iteration to the FP. 

3. |f t' a < At < 2^r' A : since D<b(AR) becomes negative, the 
orbit converges alternatingly to the FP. 

4. 2ff t' a < At: the orbit diverges. 

Therefore, any binary evolutionary code which uses an ex- 
plicit iteration scheme like Q encounters this numerical in- 
stability. This is the case even if we assume that the binary 
evolutionary code solves the equations of stellar structure with 
infinite accurracy because this instability is a result of the time- 
discretization itself which is basically unavoidable for numeri- 
cal computations. 

It is now possible to estimate how many time steps an ex- 
plicit method like Q needs at least for the computation of 
a mass transfer phase during which the maximum time step 
length Af max is used. From and (II 21 it follows that at most 
the mass fraction 



AAf ma> 
M 2 



Mi 



M 2 



At n 



- £r ^2 



(14) 



can be removed from the donor per time step. After n time steps 



the initial mass M 2 has been reduced to 



M 2 (f n ) 



1 



2 



H Y 



& - ft ^2 



M 2 (f ). 



(15) 



3 For details on maps, ve ctor fiel ds an their related nomenclature, 
the reader is referred to, e.g., Guckenheimer & Holmes 1 1983). 
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For low-mass main seque nce (MS) stars, & > 
JHiellming & WebbinkllT987l and lHiellmind[l98 
ciently small mass ratio M 2 j M\ of the b i nary, ft 
analytical approximation of IPaczvriskil J 1 97 lit , and therefore 
f j ~~ ft ~ 1- As an example: for |f = 10~ 4 which is a typical 
value for low-mass MS stars and the low-mass limits £j = -| 
and ft = — | at least 4600 time steps are necessary to reduce 
the donor mass by a factor of two according to (1 1 51 . This is the 
most optimistic limit where Af = Af max . For thermally unsta- 
ble systems, which can even approach the onset of d ynamical 
inst ability JHiellming & Webbinklll987l ISchenker etai]l2002l 
and beer & Podsiadlowski 2002), easily 10 4 or 10 5 time steps 
are necessary to reduce the donor mass by a factor of two. 

The main reason for this increased computational demand 
in the case of thermally unstable mass transfer can be under- 
stood as follows: 

The donor star of a binary in which mass transfer is ther- 
mally unstable has a deep radiative envelope. Unperturbed stars 
with a deep radiative envelope can be described (to some ex- 
tent) by a polytropic stellar structure with a polytropic index 
n < 3 and an adiabat ic index y - 5/3 (monoatomic ideal 
gas), see e.g. lHiellming & Webbink (1987), Table 1. Because 
for poly tropes of index n, = (1 — ra)/(3 - n), for radiative 
stars (with n < 3) £ s is a very large, positive number. However, 
this holds only for unperturbed stars in thermal equilibrium 
and only in the limit of inf initesimally small mass loss (see 
iHiellming & Webbinkl Jl987l) . last paragraph of Sect. II). On 
the other hand, for finite mass loss £, decreases rapidly with 
the amount of mass lost. Yet for such stars, at least during the 
initial phases of thermal timescale mass transfer, £, is still much 
larg er than unity, i.e. ty pically of order 10 - 100 (as shown in 
IHiellming & Webbinkl i 19871) . their Figs. 3 and 4). As a conse- 
quence, £j - £r is then also of the order of 10 - 100, and AM max 
as given in Eq. (114-1 is smaller and the minimum number of 
required time steps increases by that factor. 

As can be seen from dl5> . an artificial increase of Hp can 
significantly reduce the required number of time steps. This 
has been used in th e past by numerous authors to speed up the 
computations (e.g., Ham eurvlll99ll) . This approach yields the 
correct mass transfer rate only when mass transfer is close to 
stationary. However, it cannot be used if the turn-on and turn- 
off of mass transfer is important, as is the case for irradiation- 
induced mass transfer cycles (cf. Paper I). 

4.2. The Feigenbaum scenario 

The linear stability analysis discussed in Sect. 14.11 describes 
only the system dynamics near the FP, i.e., the local dynam- 
ics; we will now briefly discuss the global dynamics. 
When going to dimensionless quantities 



AR-AR 



x :— 



Hp 



_ R 2 At 



and using (0 and (0, Eq. dlOt is equivalent to 
Xn+i = x n + [1 - exp(x n )] 5 =: x n + G(x n ) =: fs(x n ). 



(16) 



(17) 
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The family of functions f$ is topologically conjugated to a fam- 
ily of 5-unimodal functions 4 : Thus, f$ is qualitatively similar 
to the well-known logistic map: 



(x) — 6x(l — x) 



(18) 



and undergoes a Feigenbaum scenario for increasing 5 which 
is called the control or chaos parameter 5 . 

For (5 = 2, corresponding to At = A? max , the FP becomes 
unstable and bifurcates into an unstable FP and a stable cycle 
of period 2. For increasing values of 6, the 2-cycle also be- 
comes unstable, and a stable 4-cycle appears, and so on. At a 
critical value of 5 » 2.7, the system dynamics become chaotic, 
and beyond this value the chaotic dynamics are permanently 
interrupted by finite windows of regular dynamics which corre- 
spond to the existence of a stable cycle of any period. Figure^ 
shows 50 iterations of dl7l of one single orbit for 5000 differ- 
ent values of 6. The FP x = corresponds to AR. At about 
6 ~ 3.1, a stable orbit of period 3 appears, whose existence 
is a mathematical proof for the existence of chaotic dynamics 
in the sys tem according to t he famous "Period three implies 
chaos" bv lLi & Yorkl ( 1975). x is a logarithm of M2 and 5 is 
linear in Af. Thus, the transitional region from stable to chaotic 
dynamics is rather small. 




Figure 1. Feigenbaum diagram of map dl7> . 50 iterations of 
one orbit for 5000 different values of 6 are shown. 

Since chaotic or "random" values in the mass transfer rate 
are not desirable, how can the chaotic dynamics be suppressed? 
As mentioned above, an artificial increase of H P will increase 
Af max and therefore shift the onset of chaos to higher values of 
5. In general, the choice of a mass transfer prescription which 



4 Unimodal functions and the ir syste m dynamics are discussed 
in detail by ICollet & Eckmanr] ^980). For an overview about 
topological conjugac y and the Feigenb a um s cenario see, e.g., 
iThompson & Stewart! Il986h and Ijacksonl ll989h . For some back- 
ground about topological conj ugacy and topo logical equivalence of 
maps, see I Arnold ll983l) and lWigginsI 4l990|). A more in-depth dis- 
cussion of map I17i can also be found in iBunindliool . 

5 fg is not topologically conjugated to a full family of <S-unimodal 
functions. Hence, unlike the logistic map, f formally does not un- 
dergo a complete Feigenbaum scenario. 



differs from (0 will only shift the onset of chaos to different 
values of 6 but it cannot suppress it. At least, if M2 grows as 
x a , fg basically keeps its behaviour (without proof) for a > 2 
and this is the case for all physical models for computing mass 
transfer. 

4.3. Controlling Chaos 

It is possible to suppress chaotic dynamics and to stabilize 
a given FP by introducing small pertubations in terms of 
"Controlling Chaos" 6 . There are two distinct approaches: the 
first one varies one system parameter by a smal l amount e„ fo r 
each time step to enforce stability of the FP llOttetalJ ll990i. 
but this method requires the a priori knowledge of the FP, and 
this is not the case for our system. The second method adds a 
small correction s n to the new iteration value, i.e., 



X n +l '•— f(Xn) + S n . 



(19) 



The simplest approa ch use s s n = K(x n - jc„_i) with a suitable 
constant K JPvragasl|l992l) . This so-called delayed dynamical 
feedback involves x n _i and x n , i.e., it uses information from 
previous iteration steps. This can extend the stability limit to 
larger values of 6 depending on K, but as numerical experi- 
ence shows, this is by far not suffici ent in our case. Th e in- 
clusion of all prior iteration values JSocolaretalJll994l) like 
s n - K(x n - jc n _i ) + /te n _i with a suitable constant R achieves 
significantly better results. For R — > 1 the stability of the FP 
can be extended to arbitrarily large values of 6 but at the ex- 
pense of an arbitrarily small basin of attraction which makes 
this approch difficult for practical application. 
A nonlinear approach of the form 



■Xn+l = f(Xn-l) + K(f(x n ) - /(X n -l)) + Rs n 



(20) 



as has been proposed bv de Sousa Vieira & Lichtenberg l(ll996) 
yields a significant improvement of the basin of attraction. For 
the specific choice of R = K and K = (Df(x) - iy l DF(x), the 
FP becomes superstable, i.e., the FP iteration i2Q\ converges 
quadratically in a sufficiently small neighbourhood of the FP. 

Since our map d!7l > is of the form/ 5 (x) = x+Gs(x), Eq. ( I20> 
is equivalent to 



x a +i = x n + DG s l (x)G$(x n ). 



(21) 



Because the position of the FP is a priori unknown, the 
best available approximation to DGg(x) is given by DGs(x n ). 
Then i2l\ turns into a standard Newton's method for Gs(x). 
Therefore, we conclude that Newton's method is most likely 

the best available method to compute x, and M2 with a reason- 
able computational effort. 



6 Th e term "Controlli ng Chaos" comes from an article of the same 
title bv lOtt et alj 1 19901) who suggested to control chaotic motion in 
nonlinear systems by small perturbations. The keyword "Controlling 
Chaos" or "Control of Chaos" later became common in that spe- 
cial field of research. In their introduction to The Control of Chaos: 
Theory and App/;'caft'o«5. lBoccaletti et allfeOOOh gave a short histori- 
cal overview of that topic. 
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5. The implicit algorithm 

While the explicit iteration scheme is a simple time integration 
scheme which requires one iteration per time step, the proposed 
implicit iteration scheme is a fixed point search which requires 
several iterations per time step until the fixed point is found 
with sufficiently high accuracy. 

Fortunately, the equations of stellar structure are typ- 
ically solved by Newton's metho d, more exp licitly by 
the so-ca lled Henyey method ( Kippenhahn et al. 1967] and 
Kip penhahn & Weigertlll990l) . Thus, the best solution is to in- 
clude the mass transfer rate or the stellar mass itself as an addi- 
tional variable in the Henyey method and to perform the fixed 
point search simultaneously with the solution of the equations 
of stellar structure (Paper I). Th i s has a lso been done indepen- 
dently bv lBenvenuto & d e Vito| d2003l) . For details of our im- 
plementation, see Biinina (2003). 

Figure|2]shows the mass transfer rate obtained with our bi- 
nary evolutionary code as a function of time for one specific 
binary system. Since our evolutionary code has not been de- 
signed for the analysis of chaotic dynamics, we have kept Af 
constant and instead have used the fact that the timescale of 
thermal relaxation varies, which is the dominant term in r' d for 
the binary system in question. When mass transfer starts, the 
thermal relaxation of the donor increases (i.e., t', decreases, c.f. 
Eq. ( I16l l). and, therefore, 6 also increases. At some point, ther- 
mal relaxation reaches a maximum and decreases afterwards, 
and so does 6. 




8.54-10 8 8.56-10 8 S.58-10 8 8.60-10 8 8.62-10 8 8.64-10 8 8.66-10 1 



Figure 2. Mass transfer rate M2[M /yr] as a function of time 
f [yr] for an early case A mass transfer in a binary system with 
Mi = M2 = 2M and using a constant time step length of 
5000 yr. The variation of the chaos parameter 6 is caused by 
the variation of thermal relaxation. The system dynamics are 
most unstable when T d is minimal. Every dot corresponds to 
one single time step in the calculation with the explicit scheme 
Q; the dashed line shows the result of the same computation, 
but using the implicit method. 

When using the explicit algorithm the resulting mass 
transfer rate undergoes a period doubling bifurcation before it 
exhibits chaotic dynamics. Within the chaotic region a peri- 
odic window of period 5 appears. After thermal relaxation and 



hence also 5 has reached its maximum, the Feigenbaum sce- 
nario evolves backwards. Since the decrease of 6 occurs slower 
than its increase before, even stable cycles of period 8 and 4 
appear before the secular mass transfer rate finally becomes 
stable. In contrast, the result obtained with our implicit algo- 
rithm for the same system and the same system parameters is 
shown by the dashed line. 

Although the implicit algorithm stabilizes the FP and even 
provides a quadratic convergence of the iteration, chaotic dy- 
namics are still present outside of a neighbourhood of the FP. 
Therefore, a good initial value for the iteration is necessary. 
We used a linear extrapolation of Mo and Hp to determine the 
initial value for M2 by (0. First, we perform typically 2-4 itera- 
tions and keep the preestimated value for M2 constant until the 
solution for the other stellar parameters has almost converged. 
Otherwise, even small fluctuations of R2 might push the next 
iteration value of M2 out of the basin of attraction of the FP 
and prevent convergence. Then, we finish with mostly 2-6 iter- 
ations using the full implicit algorithm to determine the correct 
mass transfer rate. 

Furthermore, to calculate the mass transfer rate with a nu- 
merical accuracy of better than 1%, the stellar radius has to be 
determined with a very high numerical accurracy. According 
to (0, for a scale height of Hp « 10 -4 ^2, which is typical for 
low-mass MS stars, R2 has to be determined with a relative ac- 
curacy of the order of 10~ 6 in order to get M2 with a relative 
accuracy of the order of 10~ 2 . To compute the radius with such 
a high accuracy the stellar physics and especially the equation 
of state must be a very smooth function of its variables. Every 
small discontinuity, especially in the outer layers of the star, can 
cause small jumps in the stellar radius which appear, magnified 
by the factor R2/HP, as significant jumps in the mass transfer 
rate. 



6. Summary and conclusions 

We have used a simple analytical model, aid autonomous or- 
dinary differential equation to describe the time evolution of 
the mass transfer rate. We have shown that, while the FP in the 
time-continuous system, i.e., the "physical" mass transfer rate 
is stable, the FP in the time-discretized system, i.e., the "nu- 
merical" mass transfer rate becomes unstable if the length of 
the time step Af exceeds a critical value Af max given by H2\ . 
We have estimated that even in the ideal case where Af = Af max 
at least several thousand time steps are necessary to reduce the 
donor mass in a low-mass binary system by a factor of two. For 
systems with thermally unstable mass transfer, it is even worse. 

We outline a mathematical proof that the iteration equation 
for the time-discretized system shows a behaviour similar to 
that of the logistic map and, for Af > Af max , undergoes a series 
of period doublings which finally leads to chaotic dynamics, 
i.e., to apparently random values of the computed mass transfer 
rate. The choice of a different explicit prescription to calculate 
the mass transfer rate results only in a shift of the critical time 
step length Af max which, in turn, depends on the characteristic 
scale length H at the FP AR, i.e., at the secular mass transfer 
rate. 
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In terms of "Controlling Chaos" we have briefly discussed 
several methods to stabilize the FP. Various modified itera- 
tion schemes which are equivalent to different explicit itera- 
tion schemes have been discussed in the literature, but they 
all do not show sufficient stabilization. Therefore, we suggest 
that using a different explicit iteration scheme may shift Af max 
to higher values but will not solve the problem. An implicit 
scheme, Newton's method, is the most promising solution be- 
cause it stabilizes the FP formally for At — > oo, although for 
this to be the case a sufficiently good initial value for the itera- 
tion is required. 

In practice, the implicit algorithm reduces the number of re- 
quired time steps by at least a factor of 10. Another advantage 
is that the implicit algorithm either yields the "correct" mass 
transfer rate or does not converge at all, whereas the explicit 
algorithm provides a result in every case, even if it is random. 
Therefore, in our binary evolutionary calculations we reject re- 
sults of the last time step if convergence is not reached and 
restart it with a smaller Af. 

Acknowledgements. We thank A. Weiss and H. Schlattl for providing 
their stellar evolutionary code, and H. Schlattl for valuable support 
and helpful discussions about numerics. We also thank U. Kolb for 
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